Determining constitutive models for fracture prediction in ductile materials undergoing plastic deformation

ABSTRACT

Described are methods, systems, and apparatus, including computer program products for determining material constitutive relationships of fracture in a ductile material undergoing plastic deformation. A material pressure function is determined that describes dependence of a damage rule of at least a portion of the ductile material on pressure. A material azimuthal angle function is determined that describes dependence of the damage rule of the portion of the ductile material on principle stress. The damage rule of the portion of the ductile material is determined based on the material pressure function and the material azimuthal angle function. A weakening function of the portion of the ductile material is determined based on the damage rule. Material constitutive relationships of fracture in the portion of the ductile material are determined based on the weakening function and a matrix property of the ductile material.

GOVERNMENT SUPPORT

The present invention was supported by Grant No. 123163-01 from theUnited States Navy—Office of Naval Research. The Government has certainrights in the invention.

FIELD OF THE INVENTION

The present invention relates generally to computer-based methods andapparatuses, including computer program products, for predictingfracture in a ductile material undergoing plastic deformation.

BACKGROUND

The fracture of ductile metals is an important design consideration inmany industrial applications. If excessive plastic deformation and/orfracture occur in machine parts or structural members, a system may failto perform its intended functions and worse yet, catastrophic failure ofa part or system may occur. Therefore, modes of fracture that lead tofailure must be understood in order to establish suitable failurecriteria in the design of systems.

Various models have been proposed to predict the initiation of fracture.The simplest theory, maximum-stress theory, states that failure occursat a point in the structure when the maximum principal stress σ₁ reachesa critical stress in simple tension, such as the yield stress or theultimate stress. The principle stress components are ordered such thatσ₁≧σ₂≧σ₃. The failure criterion for yielding, for example, can be statedto occur when σ₁ reaches σ_(ys) (σ₁=σ_(ys)), the yield stress. Undermaximum-stress theory, the onset of fracture is not affected by thelevel of the other principal stresses. Maximum-shear and octahedralshear stress theories, on the other hand, take into account the otherprincipal stresses. Maximum-shear stress theory, or Tresca's theory,assumes that failure occurs in a body when the maximum shear stress at apoint reaches the same value of stress for failure in a simple tensiontest. For example, maximum-shear stress theory states that yieldingfailure occurs when σ_(max)−σ_(min)=σ_(ys). Octahedral shear stresstheory, or von Mises' theory, states that failure at a particularlocation occurs when the energy of distortion reaches a critical value,such as the same energy for failure in simple tension. The failurecriterion for yielding, for example, under the von Mises theory can bestated to occur when

(σ₁−σ₂)²+(σ₁−σ₃)²+(σ₂−σ₃)²=2σ_(ys) ².  (1)

These simple threshold conditions work well with brittle materials thatfail in the linear elastic range. For ductile metals, however,irreversible plastic deformation starts when yielding occurs. Therefore,more complicated criteria are required to better model the underlyingphysics of fracture in ductile materials.

The principal stress system can be equally represented in the Cartesiancoordinate system (σ₁, σ₂, σ₃) or a cylindrical coordinate system (p, θ,σ), where p is the hydrostatic pressure, θ is the azimuthal angle on anoctahedral plane, and σ is the equivalent stress. For isotropicmaterials, the azimuthal angle can be represented by the Lode angleθ_(L), where −π/6≦θ_(L)≦π/6, because of permutation symmetry.

Both macroscopic and microscopic models have been proposed fromdifferent perspectives, such as fracture strain locus, energyabsorption, and void growth, to predict failure modes of ductilematerials. Macroscopic models generally assume homogeneous materials andare constructed based on macroscopic field variables. Unlike macroscopichomogeneous assumptions, microscopic models treat the materials asclusters of inhomogeneous cells. Microscopic models of metals are, ingeneral, more complex and can include descriptions of, for example, themultiphase materials, such as grains, precipitates, voids, and theirevolution.

One representative macroscopic fracture model is the Johnson-Cook model(G. R. Johnson, W. H. Cook, “Fracture characteristics of three metalssubjected to various strains, strain rates, temperatures and pressures,”Engineering Fracture Mechanics, 21(1):31-48, 1985). The Johnson-Cookstrength model is an empirical equation used in simulations for itsconsideration of plastic hardening, strain rate sensitivity andtemperature softening. In the Johnson-Cook model, the equivalent stress,σ, is expressed as the following function of the equivalent plasticstrain, ε_(p), and the temperature, T, and the plastic strain rate, {dotover (ε)}_(p):

$\begin{matrix}{{\overset{\_}{\sigma} = {{\left\lbrack {A + {B\; ɛ_{p}^{n}}} \right\rbrack \left\lbrack {1 + \; {C\; {\ln \left( \frac{{\overset{.}{ɛ}}_{p}}{{\overset{.}{ɛ}}_{0}} \right)}}} \right\rbrack}\left\lbrack {1 - \left( \frac{T - T_{room}}{T_{melt} - T_{room}} \right)^{q}} \right\rbrack}},} & (2)\end{matrix}$

where ε ₀ is a reference plastic strain rate, T_(room) and T_(melt) arethe room temperature and the material melting temperature, respectively,and A, B, C, n, and q are five material constants. The Johnson-Cookmodel accounts for isotropic strain hardening, strain rate hardening,and temperature softening in the uncoupled form. The first term of theright hand side in Eq. (2) represents the quasi-static stress-strainrelation at room temperature; the second term represents the strain-ratehardening; the third term represents the temperature dependence of thestress-strain relation. The fracture criterion used in the Johnson-Cookmodel is a weighted integral with respect to the effective strain:

$\begin{matrix}{{D = {\int_{0}^{ɛ_{c}}\frac{ɛ_{p}}{ɛ_{f}\left( {\frac{\sigma_{m}}{\overset{\_}{\sigma}},{\overset{.}{ɛ}}_{p},T} \right)}}},} & (3)\end{matrix}$

where

$ɛ_{f}\left( {\frac{\sigma_{m}}{\overset{\_}{\sigma}},{\overset{.}{ɛ}}_{p},T} \right)$

describes the fracture envelope, σ_(m) is the mean stress, and ε_(c) isthe plastic strain when fracture occurs.

The functional form of the weighting function is the reciprocal of theeffective failure strain, which is treated as a function of the stresstriaxiality status, the strain rate and the temperature:

$\begin{matrix}{ɛ_{f}{\quad{\left( {\frac{\sigma_{m}}{\overset{\_}{\sigma}},{\overset{.}{ɛ}}_{p},T} \right) = {\left\lbrack {D_{1} + {D_{2}{\exp \left( {D_{3}\frac{\sigma_{m}}{\overset{\_}{\sigma}}} \right)}}} \right\rbrack \left\lbrack {1 + {D_{4}{\ln \left( \frac{{\overset{.}{ɛ}}_{p}}{{\overset{.}{ɛ}}_{0}} \right)}}} \right\rbrack}}\quad}{\quad{\left\lbrack {1 - {D_{5}\; \frac{T - T_{room}}{T_{melt} - T_{room}}}} \right\rbrack.}}} & (4)\end{matrix}$

where D₁, D₂, D₃, D₄, and D₅ are material constants. The Johnson-Cookmodel does not take into account any dependence of the fracturecriterion on the Lode angle (or the third stress invariant).

Another representative macroscopic model is the Wilkins model (M. L.Wilkins, et al., “Cumulative Strain-Damage Model of Ductile Fracture,”Technical Report UCRL-53058, Lawrence Livermore National Laboratory,October 1980). Wilkins et al. proposed the failure strain envelop as aproduct of two parts: the hydrostatic tension and the stress asymmetry.The original form of Wilkins model is

$\begin{matrix}{{D = {{\int_{0}^{ɛ_{c}}{w_{1}w_{2}{ɛ_{p}}}} = 1}},} & (5)\end{matrix}$

where w₁=1/(1+ap)^(α), w₂=(2−A)^(β), A=max{s₂/s₁,s₂/s₃}, and s₁, s₂, ands₃ are the ordered principal deviatoric stress components. The materialconstants are a, α, and β. The pressure dependence is taken into accountin w₁, while w₂ accounts for the Lode angle dependence. The Wilkinsdamage model is linear with respect to the plastic strain. The Wilkinsmodel does not use a weakening or softening model in its fracturecriterion.

A representative microscopic model was proposed by Lemaitre in “AContinuous Damage Mechanics Model for Ductile Fracture,” Journal ofEngineering Materials and Technology, 107:83-89, 1985. Lemaitre'sproposed model differs from void nucleation growth coalescence (VNGC)models in that the micro-mechanism of the growth of individual voids andtheir interactions are depicted in a phenomenological way. Macroscopicparameters and equations are used in Lemaitre's model to describe theaggregative responses of the body of solids. Therefore, the constitutiveand damage model of the material is directly based on the externallyobserved behavior of the material and no microscopic interpretation ofthe material structure is required. Continuum damage mechanics considersthe material property as the combination of the matrix material anddamage accumulation. Lemaitre determined damage accumulationexperimentally for a light alloy. Their results show an increasingtendency of damage accumulation. The law of elasticity coupled withdamage is σ_(ij)=C_(ijkl)ε_(kl) ^(e)(1−D), where C_(ijkl) is the fourthorder elasticity stiffness tensor. Here, the damage variable D has ameaning of material deterioration of the stiffness. Lemaitre's modeldoes not take into account any dependence of the fracture criterion onthe Lode angle.

Another representative microscopic model is theGurson-Tvergaard-Needleman model, the so-called GTN model, which takesinto account the details of void-nucleation-growth and coalescence.Gurson derived the pressure dependent yield function for voidscontaining cell structure. See, A. L. Gurson, “Continuum theory ofductile rupture by void nucleation and growth: Part I Yield criteria andflow rules for porous ductile media,” Journal of Engineering Materialsand Technology, 99:2-15, 1977. Tvergaard and Needlemen further extendedthis model to include the void nucleation and coalescence process. TheGTN model has been used to predict several important fracture modes,such as the cup-cone fracture of a tensile round bar. See, V. Tvergaard,“Influence of voids on shear band instabilities under plane strainconditions,” International Journal of Fracture, 17:389-407, 1981, and V.Tvergaard and A. Needleman, “Analysis of the cup-cone fracture in around tensile bar,” Acta Metallurgica, 32(1):157-169, 1984. The GTNmodel does not take into account any dependence of the fracturecriterion on the Lode angle.

SUMMARY OF THE INVENTION

One approach to predicting fracture in a ductile material is todetermine material constitutive relationships of fracture in the ductilematerial undergoing plastic deformation. In one aspect, there is amethod of determining material constitutive relationships of fracture inthe ductile material undergoing plastic deformation. The method includesdetermining a material pressure function describing dependence of damageaccumulation of at least a portion of the ductile material on pressure.The method also includes determining a material azimuthal angle functiondescribing dependence of damage accumulation of the portion of theductile material on principle stress. The method further includesdetermining a damage rule of the portion of the ductile material basedon the material pressure function and the material azimuthal anglefunction. The method also includes determining a weakening function ofthe portion of the ductile material based on the damage rule. The methodalso includes determining material constitutive relationships offracture in the portion of the ductile material based on the weakeningfunction and a matrix property of the ductile material.

In another aspect, there is a method of determining materialconstitutive relationships of fracture in the ductile materialundergoing plastic deformation. The method includes determining amaterial pressure function describing dependence of damage accumulationof at least a portion of the ductile material on pressure. The methodalso includes determining a material azimuthal angle function describingdependence of damage accumulation of the portion of the ductile materialon principle stress. The method further includes determining a damagerule of the portion of the ductile material based on the materialpressure function and the material azimuthal angle function. The methodalso includes determining a weakening function of the portion of theductile material based on the damage rule. The weakening function is afunction of damage accumulation and the rate of change of the weakeningfunction with respect to damage accumulation is dependent on damageaccumulation. The method also includes determining material constitutiverelationships of fracture in the portion of the ductile material basedon the weakening function and a matrix property of the ductile material.

In another aspect, there is a computer program product. The computerprogram product is tangibly embodied in an information carrier, thecomputer program product including instructions being operable to causea data processing apparatus to determine a material pressure functiondescribing dependence of damage accumulation of at least a portion ofthe ductile material on pressure. The computer program product alsoincluding instructions being operable to determine a material azimuthalangle function describing dependence of damage accumulation of theportion of the ductile material on principle stress. The computerprogram product also including instructions being operable to determinea damage rule of the portion of the ductile material based on thematerial pressure function and the material azimuthal angle function.The computer program product also including instructions being operableto determine a weakening function of the portion of the ductile materialbased on the damage rule. The computer program product also includinginstructions being operable to determine material constitutiverelationships of fracture in the portion of the ductile material basedon the weakening function and a matrix property of the ductile material.

Another approach to predicting fracture in a ductile material is todetermine material constitutive relationships of fracture in the ductilematerial undergoing plastic deformation. In one aspect, there is amethod of determining material constitutive relationships of fracture inthe ductile material undergoing plastic deformation. The method includesdetermining a material pressure function describing dependence of adamage rule of at least a portion of the ductile material on pressure.The method also includes determining a material azimuthal angle functiondescribing dependence of the damage rule of the portion of the ductilematerial on principle stress. The method also includes determining anonlinear damage accumulation rule for the portion of the ductilematerial based on the material pressure function and the materialazimuthal angle function. The method can include determining a weakeningfunction of the portion of the ductile material based on the nonlineardamage accumulation rule. The method can also include determiningmaterial constitutive relationships of fracture in the portion of theductile material based on the weakening function and a matrix propertyof the ductile material.

Another approach to predicting fracture in a ductile material is todetermine material constitutive relationships of fracture in the ductilematerial undergoing plastic deformation. In one aspect, there is amethod of determining damage accumulation characteristics of materialconstitutive relationships of fracture in a ductile material undergoingplastic deformation. The method includes determining a material pressurefunction describing dependence of damage accumulation of at least aportion of the ductile material on pressure. The method also includesdetermining a material azimuthal angle function describing dependence ofdamage accumulation of the portion of the ductile material on principlestress. The method also includes determining a damage rule of theportion of the ductile material based on the material pressure functionand the material azimuthal angle function. The method also includesdetermining a damage accumulation in an integral form based on thematerial damage rule, the pressure function, and the material azimuthalangle function.

In other examples or embodiments, any of the aspects above can includeone or more of the following features. The method can also includedetermining an initiation of fracture in the portion of the ductilematerial based on the weakening function and a matrix property of theductile material. The material pressure function can include alogarithmic function of pressure. The material pressure function can beμ_(p)(p)=1−q log(1−p/p_(lim)). When a relative ratio of principalstresses is between about 0 and about 0.5, the material azimuthal anglefunction can be μ_(θ)=√{square root over (χ²−χ+1)}/(1+(√{square rootover (3)}/γ−2)χ), and when the relative ratio of principal deviatoricstresses is between about 0.5 and about 1, the material azimuthal anglefunction can be μ_(θ)=√{square root over (χ²−χ+1)}/(1+(√{square rootover (3)}/γ−2)(1−χ)). The material azimuthal angle function can be amaterial Lode angle function. The material Lode angle function can beμ₀=γ+(1−γ)(6|θ_(L)|/π)^(k). The material azimuthal angle function can beμ₀=1−(1−γ)(6|θ|/π)^(k) for 0≦θ≦π/6. The material azimuthal anglefunction can be μ₀=1−(1−γ)(6|π/3−θ|/π)^(k) for 0≦θ≦π/3. The damage rulecan be a function of the ratio of a plastic strain, ε_(p), to a fracturestrain, ε_(f). The damage rule can be D=(ε_(p)/ε_(f))^(m). The damagerule can be D=(exp[λ(ε_(p)/ε_(f))^(m)]−1)/(exp[λ]−1). The damage rulecan be a function of the ratio of a plastic distortion ε_(d) to afracture strain ε_(p). The damage rule can be D=(ε_(d)/ε_(p))^(m). Thedamage rule can be D=(exp[λ(ε_(d)/ε_(f))^(m)]−1)/(exp[λ]−1). Theweakening function can be a linear function of damage accumulation. Theweakening function can be w(D)=(1−D^(β))^(1/η). β can be greater thanone. η can be greater than one. An equivalent stress can be used todetermine material constitutive relationships of fracture. Theequivalent stress can be based on the weakening function and a matrixproperty of the ductile material. The equivalent stress can beσ=w(D)σ_(M).

Any of the above implementations can realize one or more of thefollowing advantages. By incorporating both the pressure dependence andthe azimuthal angle dependence into damage and weakening rules,numerical fracture prediction codes can provide simulated results withimproved accuracy with respect to the physical systems they aremodeling. By incorporating a weakening function and a damage rule,numerical simulations can maintain a continuum assumption that allowsfor a simplified model. Furthermore, numerical simulations incorporatinga weakening function and a damage rule based on pressure and azimuthalangle can provide accurate simulations of ductile materials that havegone beyond moderate plastic deformation. Under this model, fracturesimulations can be used by a wide range of industrial users (e.g.,designers of bridges, ships, vehicles, airplanes, buildings, consumerappliances, etc.) to improve the safety of products, lower developmentcosts, and/or speed up the time frame of research and design.

Other aspects and advantages of the present invention will becomeapparent from the following detailed description, taken in conjunctionwith the accompanying drawings, illustrating the principles of theinvention by way of example only.

BRIEF DESCRIPTION OF THE DRAWINGS

The foregoing and other objects, features, and advantages of the presentinvention, as well as the invention itself, will be more fullyunderstood from the following description of various embodiments, whenread together with the accompanying drawings.

FIG. 1A is a flowchart depicting the determination of materialconstitutive relationships of fracture in a ductile material undergoingplastic deformation.

FIG. 1B is a flowchart depicting the determination of materialconstitutive relationships of fracture in a ductile material undergoingplastic deformation.

FIG. 2 illustrates a load displacement curve illustrating the sequencesleading to fracture in a load carrying material.

FIG. 3 illustrates a fracture envelope represented in the threedimensional space of plastic strain and the hydrostatic tension.

FIG. 4 illustrates a hydrostatic pressure function used to fit withexperimental data.

FIG. 5 illustrates the Lode angle in the octahedral plane.

FIG. 6 illustrates a Lode angle dependence function of the first kind,which is a linear function in the plastic strain plane.

FIG. 7 illustrates a Lode angle dependence function of the second kind,which is a curvilinear function in the plastic strain plane.

FIG. 8 illustrates a fracture envelope in the principal stress spacewith the combined effects of the hydrostatic pressure and the deviatoricstate.

FIGS. 9A-9C illustrates the displacement history of a repeated loadingpath, the total equivalent plastic strain, and the equivalent plasticdistortion and the damage accumulation using equivalent plastic strainand equivalent plastic distortion for a strain controlled fractureloading.

FIG. 10 illustrates the damage rule functions to the ratio of theplastic strain to the fracture strain with damage exponent m=2 and theweakening function for several weakening exponents β and η.

FIG. 11 illustrates a hypothetical matrix stress-strain curve and theweakened material stress-strain curve for several weakening exponent βvalues with η=1.

DETAILED DESCRIPTION

The term “ductile fracture” can have multiple meanings. “Ductilefracture” can be referred to as the void nucleation-growth-coalescence(VNGC) type of fracture observed, which is usually rough. In a moregeneral sense, “ductile fracture” refers to fractures where the materialexperiences large plastic deformation while not related to a specifictype of fracture. Ductile fracture can result from the accumulation ofplastic damage.

The problem of fracture initiation in large structures involves bothmaterial and geometrical nonlinearities. Therefore, numerical approachesto fracture problems are well suited for many engineering applications.The present invention relates generally to a new method using aconstitutive material property model for ductile materials that includesfour effects: the pressure sensitivity, the azimuthal angle dependence(or, equivalently, the relative ratio of the principle stresses, and/orthe Lode angle), a damage rule, and a material weakening effect.

FIG. 1A illustrates a flowchart 100 depicting the determination ofmaterial constitutive relationships of fracture in a ductile materialundergoing plastic deformation. The determined material constitutiverelationships can be used in numerical simulations to model arbitraryductile materials for design and analysis purposes, includingdetermination of fracture criterion of parts or systems. A materialpressure function is determined (101), which can incorporate materialductility and/or fracture initiation dependence on pressure. Thispressure effect on the fracture strain can be the result of, forexample, the suppression of the initiation and propagation of microcracks and voids. The material pressure function can be determined, forexample, by fitting nonlinear equations to experimental data. A materialazimuthal angle function is determined (102), which can incorporatematerial ductility and/or fracture initiation dependence on theazimuthal angle. The azimuthal angle, which can be the Lode angle, canbe one-to-one mapped with a relative ratio of the principal stresscomponents, so that the relative ratio of the principal stresses can beused instead of the azimuthal angle. A damage rule is determined (103),which can be based on a restricted loading path, the material pressurefunction, and/or the material azimuthal angle function. The damage rulecan also be based on the ratio of the plastic strain to the fracturestrain. The damage rule, for example, can be based on a damage integralfunction. For example, incremental damage can be determined by comparingthe current plastic strain to the fracture strain at the current stressstate. A weakening function is determined (104), which can be based on adamage rule and/or damage state. The weakening function can account formaterial deterioration while maintaining a continuum assumption of thematerial. Material constitutive relationships of fracture are determined(105), wherein the strength of the matrix material can be assumed to bea basic property of the material, while the macroscopic behavior of thematerial can be affected by a weakening factor that accounts for themicroscopic damage.

FIG. 1B illustrates a flowchart 150 depicting an embodiment ofdetermining material constitutive relationships of fracture in a ductilematerial undergoing plastic deformation. A material pressure function isdetermined (151), which can incorporate material ductility and/orfracture initiation dependence on pressure. A material azimuthal anglefunction is determined (152), which can incorporate material ductilityand/or fracture initiation dependence on azimuthal angle. A damage ruleis determined (153), which can be based on a restricted loading path,the material pressure function, and/or the material azimuthal anglefunction. The damage rule can also be based on the ratio of the plasticstrain to the fracture strain. Damage accumulation is determined (161)in an integral fashion by combining the damage rule with the pressurefunction and the azimuthal angle function. A matrix strength of thematerial is determined (162), wherein the strength of the matrixmaterial can be assumed to be a basic property of the material, whilethe macroscopic behavior of the material can be affected by a weakeningfactor that accounts for the microscopic damage. A weakening function isdetermined (163), which can be based on a damage rule and/or damagestate. Mechanical properties of the material are determined (170), whichcan be based on the damage accumulation, the material matrix strength,and the material weakening function.

Plasticity Model

FIG. 2 illustrates a load displacement curve 200 that shows thesequences leading to fracture in a load carrying material. Materialfracture can be characterized, for example, by a complete loss of itsload carrying capacity or deformability. Usually, a material displayssome elasticity which is fully recoverable when the load is removed.However, beyond a certain limit, often called the elastic limit, a partof the deformation is unrecoverable. The elastic limit occurs at theyield strength 210. The region above the yield strength 210 is oftencalled plasticity. Due to strain hardening, deformation can still bestable until the ultimate strength 220 is reached. Beyond the ultimatestrength 220, the deformation is unstable and fracture 230 can soonresult. It can be assumed that the matrix stress and the plastic straincurve is a basic material property,

σ_(M)=σ_(M)(ε_(p)),  (6)

where σ_(M) is the matrix stress and ε_(p) is the plastic strain.Material deterioration can be characterized by a weakening factor, whichcan be a function of the damage, e.g.,

σ=w(D)σ_(M),  (7)

where σ is the applied stress and w(D) is a weakening or softeningfunction defined on the damage variable D.

Up to moderate plastic deformation, the microstructure of materials canbe viewed as being undamaged. The strength of material is oftensimplified as a one-dimensional problem where the resistance of materialis characterized as a function of the equivalent plastic strain only.Classical plasticity theory, or conventional continuum mechanics, isacceptable in practical problems where only small ductile damage isinvolved. However, these theories are not sufficient in dealing withmany engineering applications where deformation goes beyond moderateplasticity, when changes in microstructure may no longer beinsignificant and cannot be ignored in order to obtain accurate fractureprediction. In such conditions, the ignorance of the change ofmicrostructure is over simplified. A damage rule and a damage inducedweakening effect have to be included in the model to characterize thematerial deterioration.

From experimental observations, ductile failure of structures comprisesthree phases: (1) accumulation of damage, (2) initiation of fracture,and (3) crack propagation. Fracture initiation can be modeled as theresult of the accumulation of the ductile plastic damage.Microscopically, such damage can be associated with void nucleation,growth and coalescence, shear band movement, propagation ofmicro-cracks, and/or the like. Macroscopically, the degradation of thematerial can be exhibited as a decrease in the material stiffness,material strength, and/or a reduction of the remaining ductility. Thesephysical changes can be used as an indicator to predict the onset offracture, either based on the current value or in a cumulative fashion.

Various experiments have shown that the damage process that eventuallyleads to the fracture of a material exhibits strong dependence on theloading history of the material (e.g., the damage). Therefore, the modelincludes a damage variable that takes into account the effect of amaterial's loading history. The damage evolves when the solid materialis subjected to plastic loadings, which can be treated in the threedimensional space of the principal stresses. The principal stresses canbe represented in the cylindrical coordinate system denoted by (p, θ,σ), where p is the hydrostatic pressure, θ is the azimuthal angle on theoctahedral plane, and u is the von Mises equivalent stress. On anoctahedral plane, the pressure is fixed. The aximuthal angle can befixed on the octahedral plane, which represents a class ofdeviatorically proportional loading path. The damage evolution can thendepend only on the ratio of the plastic strain to the fracture strainfor monotonic loadings, i.e., D=f(ε_(p)/ε_(f)), where D is theaccumulated damage on the path, ε_(p) is the plastic strain and ε_(f) isthe fracture strain on that path. This path can be called the“restricted loading path” for convenience. The rate form can berepresented as

$\begin{matrix}{\overset{.}{D} = {\frac{\partial{f\left( \frac{ɛ_{p}}{ɛ_{f}} \right)}}{\partial\left( \frac{ɛ_{p}}{ɛ_{f}} \right)}{\frac{{\overset{.}{ɛ}}_{p}}{ɛ_{f}}.}}} & (8)\end{matrix}$

In Eq. (8), the normalization denominator ε_(f) is defined by the givenpressure and the azimuthal angle θ, i.e., ε_(f)=ε_(f)(p,θ). Assuming theeffects of the hydrostatic pressure and the azimuthal angle areindependent of each other, the equivalent fracture strain can be writtenas

ε_(f)=ε_(f0)μ_(p)(p)μ_(θ)(θ),  (9)

where ε_(f0) is a reference fracture strain, μ_(p)(p) represents theeffect of the hydrostatic pressure, and μ_(θ)(θ) represents thedeviatoric state. When both μ_(p)(p) and μ_(θ)(θ) are unity, the modeldegenerates to the constant failure strain criterion. The fractureenvelope defined by ε_(f) shows the relative extent of ductility for thegiven pressure and the azimuthal angle. At a more accurate level, theactual damage is an accrued effect of the loading history of varyingpressure and azimuthal angle, which can be described by a historyvariable.

Based on these assumptions, the constitutive relationship of thematerial can be described by a set of four equation, Eqs. (6)-(9). Theseequations can form a theoretical basis of a new damage plasticity model.The internal variables are the plastic strain and the damage. The inputfunctions for the material are the five curves: σ_(M)(ε_(p))=thestress-strain curve for the matrix material, w(D)=weakening function,D(ε_(p)/ε_(f))=damage rule, μ_(p)(p)=material pressure function, andμ_(θ)(θ)=material azimuthal angle function. This method to calculate thedamage is called “cylindrical decomposition.”

In this application, “cylindrical decomposition” refers to the method ofusing the pressure, the aximuthal angle (or, equivalently, the Lodeangle), and the plastic strain, which forms a cylindrical coordinatesystem in a three-dimensional space. The pressure sensitivity and theazimuthal angle dependence of fracture are used to construct a fracturesurface oriented in the hydrostatic axis in the principal stress space.The fracture surface can be equally represented in the space of plasticstrain and the hydrostatic tension, as illustrated in FIG. 3. A fracturesurface can be determined based on the hydrostatic pressure p and therelative ratio of the principal deviatoric stresses: χ=(s₂−s₃)/(s₁−s₃),which can be used interchangeably with the Lode angle. The principaldeviatoric stresses, s_(i), can be defined ass_(i)=σ_(i)−σ_(m)=σ_(i)−(σ₁+σ₂+σ₃)/3, where σ_(i) are the principalstresses and i=1, 2, 3.

Pressure Sensitivity

One approach to constructing the fracture surface includes usingexperimental data to determine pressure dependence and azimuthal angledependence on fracture. From the orthogonality of p and χ, it can beapproximated that the effects of p and χ are independent of each other.Thus, the effects of p and χ can be uncoupled to a first approximationand can be quantified separately.

Experimental data used to determine the pressure and the azimuthal angledependence on the fracture surface can be obtained by adjusting thepressure and the azimuthal angle separately, such that two series oftests can be conducted at constant p and constant χ, respectively. Forexample, some experiments show that materials can become more ductile asthey experience higher compressive pressures. For some metals, theductility can be an order higher under high compressive pressure thanits ductility at atmospheric pressure. For example, zinc exhibits asudden transition from brittle to ductile at a critical pressure ofabout 70 MPa. One test, for example, that can be used to determineexperimental data is the unidirectional tension test with confiningpressure.

Round bars can be used as specimens and pulled in a pressure chamber totest the effect of hydrostatic pressure on the material fracture strain.The round bar specimens can be viewed as being in a uniaxial tensioncondition superimposed by a hydrostatic pressure up to a compressurepressure of about, for example, 30 kbar. Because the two lateralprinciple stress components remain identical at the symmetric line forthe axisymmetric specimen, the ratio of the principal stresses remainsessentially zero at the center of the neck throughout these experiments.The ratio of the cross-sectional area at the neck at fracture to theinitial cross-sectional area usually decreases with respect to thelateral confining pressure.

A series of experimental tests have been performed on tensile round barson armor steels using a pressure chamber (Bridgman, P. W., Studies inlarge plastic flow and fracture. McGraw-Hill Inc., 1952.). From the testresults, the relationship of the fracture strain and the confiningpressure can be expressed as

$\begin{matrix}{{\frac{A_{f}}{A_{0}} = {\frac{A_{f\; 0}}{A_{0}}\left( {1 - \frac{p_{conf} - p_{atm}}{p_{\lim} - p_{atm}}} \right)^{\overset{\_}{q}}}},} & (10)\end{matrix}$

where q is a material constant that fits experimental data, A₀ is theoriginal cross-sectional area, A_(f) is the cross-sectional area afterfracture at the neck, A_(f0) is the cross-sectional area after fractureat the neck at atmospheric pressure, p_(atm) is the atmosphericpressure, and p_(lim) is a limiting pressure beyond which the materialwill not fail in the uniaxial tensile condition. The limiting pressurehas been estimated, for example, to be between about 1 GPa to about 2GPa for some steels.

Assuming p_(atm)<<p_(lim) and taking the logarithm on both sides of Eq.(10), the following equation can be obtained

$\begin{matrix}{{ɛ_{f} = {ɛ_{f\; 0}\left\lbrack {1 - {q\; {\log \left( {1 - \frac{p}{p_{\lim}}} \right)}}} \right\rbrack}},} & (11)\end{matrix}$

where ε_(f) is the fracture strain at the confining pressure p, ε_(f0)is a reference failure strain at zero mean stress, ε_(f0)=log(A₀/A_(fo))is the uniaxial tensile failure strain without confining pressure, q=isa shape parameter, and p_(lim) is a limiting pressure beyond which nodamage occurs. Therefore, the pressure dependence function μ_(p)(p)becomes

$\begin{matrix}{{\mu_{p}(p)} = {1 - {q\; {{\log \left( {1 - \frac{p}{p_{\lim}}} \right)}.}}}} & (12)\end{matrix}$

For numerical implementation, the confining pressure can be replaced bythe hydrostatic pressure and the form of Eq. (12) is retained, exceptthat the material constants can be re-calibrated for the hydrostaticpressure.

The average hydrostatic pressure 410 experienced in the course of apulling test can be illustrated as in FIG. 4. The pulling test starts at420 with no deformation. The equivalent stress path on the σ−p plane isshown as the thick solid line 430. Line 440 is the fracture loci of σwith respect to the hydrostatic pressure. The mean hydrostatic pressurealong the entire loading path can be estimated to bep_(ave)=p_(conf)−σ_(flow)/6. As illustrated in FIG. 4, the pressure atthe fracture point is p_(conf)−σ_(flow)/3.

Azimuthal Angle Dependence

From the pressure point of view, the simple shear is zero mean stressand the simple tension is positive mean stress. Experiments have shownthat simple shear (or torsion) fracture strain can be less than thesimple tension fracture strain for at least some materials. Thisexperimental result can not be explained by the pressure dependence.Therefore, an azimuthal dependence function characterizing the azimuthalvariation of the fracture locus on the same pressure can be necessary inmodeling ductile fracture.

For isotropic materials, the principal stresses are interchangeable toreflect the independence of damage to the observation frame. Theazimuthal angle dependence function for damage evolution can havepermutation symmetry in all three principal planes. Therefore, on anoctahedral plane, the azimuth angle can be divided into six parts thathave the same weighting function. Each part covers the complete range ofthe relative ratio of the principal stresses from 0 to 1. Assuming thatthe forward motion and the backward motion induce the same amount ofdamage, an additional symmetry can be introduced to the octahedralplane—i.e., the shape of the twelve segments are the same (apart fromthe reflections).

For isotropic materials, the azimuth angle can be characterized by theLode angle, which can be defined as

$\begin{matrix}{{\theta_{L} = {\tan^{- 1}\left\{ {\frac{1}{\sqrt{3}}\left\lbrack {{2\left( \frac{s_{2} - s_{3}}{s_{1} - s_{3}} \right)} - 1} \right\rbrack} \right\}}},} & (13)\end{matrix}$

where s₁, s₂, and s₃ are the maximum, intermediate and minimum principaldeviatoric stresses respectively. FIG. 5 illustrates the Lode angle. Fora given principal deviatoric stress represented by vector O'A, the Lodeangle is represented by θ_(L).

The relative ratio of the principal deviatoric stresses, χ, can bedefined as χ=(s₂−s₃)/(s₁−s₃). To recognize the difference in planestrain and triaxial tension, a new material parameter, γ, can be definedas the ratio of the fracture strain at the plane strain condition tothat at the triaxial tension (or compression) condition at the sameconstant hydrostatic pressure, i.e.

$\begin{matrix}{\gamma = {\frac{ɛ_{f}\left( {\chi = 0.5} \right)}{ɛ_{f}\left( {\chi = 0} \right)}.}} & (14)\end{matrix}$

FIG. 6 illustrates an example of a Lode dependence function that is afirst order linear relationship in the plane of plastic straincomponents. The fracture point representing the fracture strain at thetriaxial tension or compression is connected to that of the plane straincondition by a straight line and, thus, forms a polygon on the strainplane. Using the relative ratio of the stress deviators, the “six pointstar” illustrated in FIG. 6 can be represented by the function

$\begin{matrix}{\mu_{\theta} = \left\{ {\begin{matrix}{\frac{\sqrt{\chi^{2} - \chi + 1}}{1 + {\left( {\frac{\sqrt{3}}{\gamma} - 2} \right)\chi}},} & {{0 \leq \chi \leq 0.5};} \\{\frac{\sqrt{\chi^{2} - \chi + 1}}{1 + {\left( {\frac{\sqrt{3}}{\gamma} - 2} \right)\left( {1 - \chi} \right)}},} & {{0.5 \leq \chi \leq 1};}\end{matrix},} \right.} & (15)\end{matrix}$

which is symmetric with respect to χ=0.5. Therefore, a fracture envelopedefining a region of material properties in which fracture will occurcan be constructed for each of the twelve pie slices on the octahedralplane, which are identified by either 0≦χ≦0.5 or 0.5≦χ≦1. The first kindof the Lode dependence function reduces to a right hexagon whenγ=√{square root over (3)}/2. The axes of the right hexagon is ameasurement of strain.

Experimental data shows that for many materials, the material parameterγ which governs the difference in ductility between shear and tension atthe same pressure is less than unity. Consequently, the fracture favorsa shear mode. This is consistent with experimental data that shows theformation of a shear lip at the edge of an otherwise tensile specimen.

FIG. 7 illustrates another example of an azimuthal angle dependencefunction, which is a curvilinear model defined by

$\begin{matrix}{\mu_{\theta} = \left\{ {\begin{matrix}{{1 - {\left( {1 - \gamma} \right)\left( \frac{6{\theta }}{\pi} \right)^{k}}},} & {{0 \leq \theta \leq \frac{\pi}{6}};} \\{{1 - {\left( {1 - \gamma} \right)\left( \frac{6{{{\pi/3} - \theta}}}{\pi} \right)^{k}}},} & {{\frac{\pi}{6} \leq \theta \leq \frac{\pi}{3}};}\end{matrix},{{{or}\mspace{14mu} \mu_{\theta}} = {\gamma + {\left( {1 - \gamma} \right)\left( \frac{6{\theta_{L}}}{\pi} \right)^{k}}}},} \right.} & (16)\end{matrix}$

where θ is the azimuthal angle with respect to the s₁ axis, θ_(L) is theLode angle, k is the Lode dependence exponent, γ=ε_(f)(θ=π/6)/ε_(f)(θ=0)is the same material constant as in the polygon model, and β is the Lodedependence exponent, or shape parameter. When k=1 in Eq. (16), thecurvilinear Lode dependence function is an Archimedes' spiral in thepolar coordinate system, which is linear with respect to the azimuthangle on the triaxial plastic strain plane. When γ=1, the curvilinearfunction degenerates to a perfect circle. Other nonlinear forms may beobtained by fitting the experimental data.

FIG. 8 illustrates the fracture envelope in the principal stress spacerepresenting the entire family of the restricted loading paths with thecombined effects of the hydrostatic pressure and the deviatoric state.The limiting hydrostatic pressure plane 80 can correspond toexperimental data that suggests that no damage occurs when thehydrostatic pressure is sufficiently high. For example, when thehydrostatic pressure is above the hydrostatic pressure plane 80, newbonds can be created at room temperature. The cutoff pressure 81represents the point in which the equivalent fracture strain ε_(fo) iszero. As illustrated in FIG. 8, the fracture locus can shrink to asingle point 81 at the triad axis. Line 83 represents the fracture lociof equivalent stress at uniaxial tensile condition with respect tohydrostatic pressure. For example, line 85 represents a referencefailure stress on the π-plane 84 under uniaxial tensile condition.Polygon 89 illustrates the failure loci on the π-plane 84. Line 86represents another deviatoric state under uniaxial tensile condition ona different constant pressure plane 87. Polygon 88 illustrates, forexample, the failure strain at different deviatoric states on theconstant pressure plane 87.

Damage Rule

As described above, in continuum damage mechanics, materialdeterioration can be described by an internal damage variable. Damagecan be modeled as an anisotropic quantity by a tensor in general.However, damage can also be assumed to be an isotropic process andtreated as a scalar. In many applications, modeling damage as anisotropic process still provides for accurate predictions. To utilizecumulative damage as a criterion to predict the onset of fracture, therelationship of damage with respect to one or more measurable quantitiescan be established based on measurable field variables.

The physical properties of a material (e.g., the elastic modulus, theductility, the local mass density, and/or the like) can change with theaccumulation of the damage. In this sense, there can be multiple ways ofmeasuring the quantity of damage. One of these ways can be viewed fromthe point-of-view of the relative deformability of the material, whichcan be obtained by comparing the deformability of the current state ofthe material with a prior state of the material, such as the initialundamaged state. In this way, ductile damage can be defined as therelative loss of deformability of the material. The relative loss ofdeformability can be expressed formally, in one example, by D=1/N*100%,where N is the number of times that the material can survive the sameloading. For example, if a material fractures after 10 times of repeatedloading, the material loses 10% of its initial deformability after theeach loading.

A macroscopic crack can result from damage accumulation. Unlike theobservation of crack initiation, damage accumulation may not be easilytracked along a deformation path. By conventional means, a single testgives only one point of the ultimate fracture in the damage accumulationgraph with respect to the plastic strain, i.e. D=1 at ε_(p)=ε_(f).

The loss of relative deformability can also be viewed by usinglow-cycle-fatigue test results. For materials under strain reversals,the total plastic strain at fracture can be greater than that ofmonotonic loading condition. This observation suggests relatively smalldamage at low plastic distortion for the same amount of incrementalequivalent strain.

Reversed deviatoric loading can induce substantially the same amount ofdamage as forward loading. Reverse deviatoric loading means that theload is at the opposite direction of the forward loading or theazimuthal angle is at 180° from forward loading, while the pressureremains the same. For example, reverse loading for uniaxial tensionσ₁=σ₀>0, σ₂=σ₃=0 is σ₁=σ₂=2σ₀/3>0, σ₃=−σ₀/3<0. Reverse loading forsimple shear σ₁=σ₀>0, σ₂=0, σ₃=−σ₀<0 is the same as its forward loading,i.e. σ₁=σ₀>0, σ₂=0, σ₃=−σ₀<0. Here, σ₀ is a positive scalar.

Assuming that reversed deviatoric loading induces the same amount ofdamage as forward loading, the fracture locus on a hydrostatic plane hasreverse symmetry in addition to permutation symmetry. The fracturestrain for triaxial tension can be the same as the triaxial compressionfor a monotonic loading at identical constant pressure. Here, the termstriaxial tension (i.e., θ_(L)=−30°) and triaxial compression (i.e.,θ_(L)=30°) should be distinguished from tension and compression, whichare used for “simple tension” and “simple compression” sometimes in theliterature. This results in an identical twelve pieces in the azimuthaldirection on a hydrostatic plane. Any point of fracture loci (except theorigin and those on the triad axes) has eleven image points.

The repeated plastic loading path (or low-cycle-fatigue test) reversesthe direction of straining at intervals. This type of loading can becharacterized by the ratio of the minimum and the maximum strain, i.e.R=ε_(min)/ε_(max). For one type of reversed loading path, R=0. A fullcycle can be viewed as two branches: a loading branch and a reversedloading branch, between which a strain reversal occurs. The plasticstrain on the current branch can be denoted by the plastic distortionε_(d) to distinguish from the total plastic strain ε_(p).

The plastic strain is a history variable which can be defined as

${ɛ_{p} = {\sqrt{\frac{2}{3}}{\int\sqrt{\left( {ɛ_{1}} \right)^{2} + \left( {ɛ_{2}} \right)^{2} + \left( {ɛ_{3}} \right)^{2}}}}},$

where dε₁, dε₂, and dε₃ are the principal plastic strain increments. Theplastic distortion is a state variable which can be defined as

${ɛ_{d} = {\sqrt{\frac{2}{3}}\sqrt{\left( ɛ_{1} \right)^{2} + \left( ɛ_{2} \right)^{2} + \left( ɛ_{3} \right)^{2}}}},$

where ε₁, υ₂, and ε₃ are the current principal plastic straincomponents. FIG. 9B illustrates the difference between the plasticstrain and the plastic distortion.

FIGS. 9A, 9B, and 9C illustrate the time dependence of the rotationalangle, the total equivalent plastic strain and the equivalent plasticdistortion, respectively. To remove the pressure dependence, a torsiontest, where pressure is constant zero, can be used to demonstrate theevolution of damage using the plastic strain and the plastic distortion.FIGS. 9A-9C show two identical cycles. FIG. 9A illustrates a rotationalrepeated plastic loading path with R=0. A loading branch 910 changes therotational angle from zero to Δθ_(max) and a reverse loading branch 920returns back to zero. FIG. 9B shows the evolution of the plastic strain930 and the plastic distortion 940. Assuming a power law damage rulewith damage exponent m=2, the damage accumulation using the plasticstrain can be plotted as the solid line 950. The damage accumulationusing the plastic distortion can be plotted as the dash-dot line 960.For this example, the predicted total fracture strain using these twomethods can be as much as twice for repeated loading conditions, asshown in FIG. 9C.

For a low cycle fatigue test, a particular form of the damage rule canbe derived from the Manson-Coffin's relationship. Experiments show thatlow-cycle fatigue is plasticity dominated phenomenon and appears to bedependent mainly upon the ductility of metals. Using the Palmgren-Minerrule, which assumes a linear relationship of damage accumulation oflow-cycle fatigue. Therefore, damage is only related to the plasticdeformation on the current branch.

Generally speaking, the number of cycles for low-cycle fatigue can beless than 10000 cycles. The relationship between the applied plasticstrain and the number of cycles to failure (Δε_(p)−N curve) can bedescribed by the Manson-Coffin relationship for a number of materials,i.e., Δε_(p)·N^(k′)=C, where C and k′ are material constants. TheManson-Coffin relationship appears to be linear on a log-log scale plotof Δε_(p) and N. Experimental data from push-pull tests and torsiontests of numerous materials such as bars of aluminum alloy, copper,and/or carbon steel show good agreement with this linear relationship.

In each loading branch, the loading is monotonic, dΔε_(p) is the same asdε_(p). The relationship of incremental damage with respect toincremental equivalent plastic strain is

$\begin{matrix}{\frac{D}{ɛ_{p}} = {\frac{1}{2k^{\prime}C^{({1/k^{\prime}})}ɛ_{p}^{({1 - {1/k^{\prime}}})}}.}} & (17)\end{matrix}$

For N=½, the material fails at a plastic strain of ε_(f) from amonotonic test. Thus, C=ε_(f)/2^(k′).

Letting m=1/k′, Eq. (17) can be rewritten as

$\begin{matrix}{{dD} = {{m\left( \frac{ɛ_{p}}{ɛ_{f}} \right)}^{({m - 1})}\frac{1}{ɛ_{f}}d\; {ɛ_{p}.}}} & (18)\end{matrix}$

Eq. (18) can be used in an integral calculation to determine damage.

“Damage accumulation” can refer to the integral equation, for example

$\begin{matrix}{{D = {\int_{0}^{ɛ_{c}}{{m\left( \frac{ɛ_{p}}{ɛ_{f}} \right)}^{({m - 1})}{ɛ_{p}}}}},} & (19)\end{matrix}$

where ε_(c) denotes a critical plastic strain at the onset of ductilefracture. A “damage rule” can refer to the integrated function, forexample, D(ε_(p)/ε_(f))=(ε_(p)/ε_(f))_(m).

For monotonic proportional plastic loading paths, the plastic strain andthe plastic distortions are identical. For repeated loading conditions,as illustrated in FIG. 9, the plastic distortion is used in the damageaccumulation function,

$\begin{matrix}{D = {\int_{0}^{ɛ_{c}}{{m\left( \frac{ɛ_{d}}{ɛ_{f}} \right)}^{({m - 1})}{{ɛ_{d}}.}}}} & (20)\end{matrix}$

to depict the actual damage rate. The power law damage rule function isillustrated in FIG. 10 for m=2 along with the weakening function forseveral weakening exponents β. Another damage rule function is theexponential damage rule function. By more careful examination of theΔε_(p)−N curve on the log-log scale plot, the curve is concave towardthe Δε_(p) axis. This damage rule can be more accurately represented by

$\begin{matrix}{{D = \frac{{\exp \left\lbrack {\lambda \left( \frac{ɛ_{p}}{ɛ_{f}} \right)}^{m} \right\rbrack} - 1}{{\exp \lbrack\lambda\rbrack} - 1}},{{{or}\mspace{14mu} D} = {\frac{{\exp \left\lbrack {\lambda \left( \frac{ɛ_{d}}{ɛ_{f}} \right)}^{m} \right\rbrack} - 1}{{\exp \lbrack\lambda\rbrack} - 1}.}}} & (21)\end{matrix}$

where λ is a material constant.

It is postulated that the damage potential is a fundamental materialbehavior. The hydrostatic pressure and the relative ratio of the stressdeviators enter into the damage potential implicitly through theequivalent fracture strain and do not change the damage potential inother form.

Material Deterioration

In continuum damage mechanics (CDM), damage is one of the constitutivevariables that takes into account the degradation and loss of loadcarrying capacity of materials. Considering a unitary reference volumeelement, the effective load carrying area can be viewed as decreasing asdamage accumulates. The loss of the effective resisting area can be usedin continuum damage mechanics to describe the material deterioration.The loss of effective resisting area can be defined as

$\begin{matrix}{{D_{s} = {1 - \frac{A_{eff}}{A_{0}}}},} & (22)\end{matrix}$

where A₀ is the nominal sectional area of the reference element andA_(eff) is the effective resisting area. The concept of A_(eff) isuseful when considering microstructural change in the solid.

For an undamaged material, no defects can be assumed to exist and,therefore, D_(s)=0. At a fully damaged state, there is no resistancearea, i.e. A_(eff)=0, or D_(s)=1. Thus, by using the strain equivalencehypothesis, the macroscopic behavior of the solid can be calculated fromthe matrix strength from the force balance of the unitary element, i.e.

$\begin{matrix}{\overset{\_}{\sigma} = {\frac{A_{eff}}{A_{0}}{\sigma_{M}.}}} & (23)\end{matrix}$

Applying Eq. (22) to Eq. (23), the equivalent strength of themacroscopic solid is

σ=(1−D _(s))σ_(M).  (24)

where a weakening factor w(D)=1−D_(s) is applied on the matrix strength.A relationship between D and D_(s) can be established to couple thedamage accumulation and the material deterioration. FIG. 11 illustratesa hypothetical matrix stress-strain curve and the weakened materialstress-strain curve for several weakening exponent β values. The otherexponent η is assumed to be unity here, i.e., η=1, and the pressureeffect is neglected. The equivalent fracture strain for thishypothetical material is 1.

The damage induced material weakening can be introduced into theconstitutive model by coupling of the yield function and the associatedflow rule with the damage. The material deterioration can be consideredat a material volume level. Following the continuum damage mechanics,the constitutive equation of the damaged material can be derived fromthe modified yield function, viz.

Φ= σ ² −[w(D _(s))σ_(M)]².  (25)

The strain for a damaged material can be represented by constitutiveequations of the undamaged material in the yield function of which thestress is simply replaced by the effective stress. The materialhardening is described by the hardening behavior of the matrix material,which is represented by σ_(M) in Eq. (25), while the material softeningcan be represented by the multiplication factor of w(D).

The stress-strain relationship of the matrix material can be assumed tobe a basic material property, i.e., it does not change with the loadingcondition and other internal state variables, such as, for example, thedamage, the pressure and the azimuthal angle. Upon accumulation of themicro structural change along with the plastic deformation, the externalbehavior the solid that contains the matrix material and the microstructural damage also changes—often, for example, in a weakened manner.Assuming the yielding of the damaged material follows von Mises type ofyield criterion, the material weakening can refer to, for example, themacroscopic material yield strength, i.e. σ≦σ_(M).

Because the damage can be evaluated in various ways (e.g., reduction offracture strain, reduction of stiffness, increase of porosity, etc.), ingeneral, the accumulation of these different types of damage may notnecessarily be the same. For example, the damage associated with thereduction of the effective load carrying area may not necessarily be thesame as the damage that addresses the reduction of the deformability,i.e. D_(s)<D or D_(s)>D.

The undamaged state can be assumed to have D=0 and a fully damage stateof D=1. The undamaged state can have full strength of the matrixmaterial, therefore, w(D=0)=1. On the other hand, the fully damagematerial no long carries any loading, therefore, w(D=1)=0. These twoequations can be the end conditions that a weakening function has tosatisfy.

For example, a class of the weakening function can be assumed to bew(D)=(1−D^(β))^(1/η), where β and η are weakening exponents. It followsthat w(D)^(η)+D^(β)=1. It can be verified that this class of weakeningfunction satisfies the end conditions. FIG. 10 illustrates weakeningfunctions associated with the ductile damage. Other nonlinearrelationships that satisfy the end condition and link the reduction ofthe effective loading carrying area to a general quantity of damage toreplace the multiplication factor (1−D_(s)) are also possible.

When β=1 and η=1, this class of weakening function reduces to a specialcase of linear weakening function with respect to the ductile damage.The stiffness damage and the ductile damage can be approximate to be thesame a first order accuracy. i.e. D_(s)=D. Experiments show that thestiffness damage accumulate can lag behind the accumulation of ductiledamage. Since Dε[0,1] and D_(s)ε[0,1], then D_(s)≦D when β>1 and η>1.

The above-described techniques can be implemented in digital electroniccircuitry, or in computer hardware, firmware, software, or incombinations of them. The implementation can be as a computer programproduct, i.e., a computer program tangibly embodied in an informationcarrier, e.g., in a machine-readable storage device or in a propagatedsignal, for execution by, or to control the operation of, dataprocessing apparatus, e.g., a programmable processor, a computer, ormultiple computers. A computer program can be written in any form ofprogramming language, including compiled or interpreted languages, andthe computer program can be deployed in any form, including as astand-alone program or as a subroutine, element, or other unit suitablefor use in a computing environment. A computer program can be deployedto be executed on one computer or on multiple computers at one site. Theabove-described techniques can be incorporated in a numerical simulationsoftware package using a finite-element method solver. Software packagescan include, for example, LS-DYNA and ABAQUS.

Method steps can be performed by one or more programmable processorsexecuting a computer program to perform functions of the invention byoperating on input data and generating output. Method steps can also beperformed by, and an apparatus can be implemented as, special purposelogic circuitry, e.g., an FPGA (field programmable gate array) or anASIC (application-specific integrated circuit). Subroutines can refer toportions of the computer program and/or the processor/special circuitrythat implements that functionality.

Processors suitable for the execution of a computer program include, byway of example, both general and special purpose microprocessors, andany one or more processors of any kind of digital computer. Generally, aprocessor receives instructions and data from a read-only memory or arandom access memory or both. The essential elements of a computer are aprocessor for executing instructions and one or more memory devices forstoring instructions and data. Generally, a computer also includes, orbe operatively coupled to receive data from or transfer data to, orboth, one or more mass storage devices for storing data, e.g., magnetic,magneto-optical disks, or optical disks. Data transmission andinstructions can also occur over a communications network. Informationcarriers suitable for embodying computer program instructions and datainclude all forms of non-volatile memory, including by way of examplesemiconductor memory devices, e.g., EPROM, EEPROM, and flash memorydevices; magnetic disks, e.g., internal hard disks or removable disks;magneto-optical disks; and CD-ROM and DVD-ROM disks. The processor andthe memory can be supplemented by, or incorporated in special purposelogic circuitry.

To provide for interaction with a user, the above described techniquescan be implemented on a computer having a display device, e.g., a CRT(cathode ray tube) or LCD (liquid crystal display) monitor, fordisplaying information to the user and a keyboard and a pointing device,e.g., a mouse or a trackball, by which the user can provide input to thecomputer (e.g., interact with a user interface element). Other kinds ofdevices can be used to provide for interaction with a user as well; forexample, feedback provided to the user can be any form of sensoryfeedback, e.g., visual feedback, auditory feedback, or tactile feedback;and input from the user can be received in any form, including acoustic,speech, or tactile input.

The above described techniques can be implemented in a distributedcomputing system that includes a back-end component, e.g., as a dataserver, and/or a middleware component, e.g., an application server,and/or a front-end component, e.g., a client computer having a graphicaluser interface and/or a Web browser through which a user can interactwith an example implementation, or any combination of such back-end,middleware, or front-end components. The components of the system can beinterconnected by any form or medium of digital data communication,e.g., a communication network. Examples of communication networksinclude a local area network (“LAN”) and a wide area network (“WAN”),e.g., the Internet, and include both wired and wireless networks.

The computing system can include clients and servers. A client and aserver are generally remote from each other and typically interactthrough a communication network. The relationship of client and serverarises by virtue of computer programs running on the respectivecomputers and having a client-server relationship to each other.

One skilled in the art will realize the invention may be embodied inother specific forms without departing from the spirit or essentialcharacteristics thereof. The foregoing embodiments are therefore to beconsidered in all respects illustrative rather than limiting of theinvention described herein. Scope of the invention is thus indicated bythe appended claims, rather than by the foregoing description, and allchanges that come within the meaning and range of equivalency of theclaims are therefore intended to be embraced therein.

1. A method of determining material constitutive relationships offracture in a ductile material undergoing plastic deformation, themethod comprising: determining a material pressure function describingdependence of damage accumulation of at least a portion of the ductilematerial on pressure; determining a material azimuthal angle functiondescribing dependence of damage accumulation of the portion of theductile material on principle stress; determining a damage rule of theportion of the ductile material based on the material pressure functionand the material azimuthal angle function; determining a weakeningfunction of the portion of the ductile material based on the damagerule; and determining material constitutive relationships of fracture inthe portion of the ductile material based on the weakening function anda matrix property of the ductile material.
 2. The method of claim 1further comprising determining an initiation of fracture in the portionof the ductile material based on the weakening function and a matrixproperty of the ductile material.
 3. The method of claim 1 wherein thematerial pressure function includes a logarithmic function of pressure.4. The method of claim 3 wherein the material pressure function is:${\mu_{p}(p)} = {1 - {q\; {{\log \left( {1 - \frac{p}{p_{\lim}}} \right)}.}}}$5. The method of claim 1 wherein when a relative ratio of principalstresses is between about 0 and about 0.5, the material azimuthal anglefunction is:${\mu_{\theta} = \frac{\sqrt{\chi^{2} - \chi + 1}}{1 + {\left( {\frac{\sqrt{3}}{\gamma} - 2} \right)\chi}}};$and when the relative ratio of principal deviatoric stresses is betweenabout 0.5 and about 1, the material azimuthal angle function is:$\mu_{\theta} = {\frac{\sqrt{\chi^{2} - \chi + 1}}{1 + {\left( {\frac{\sqrt{3}}{\gamma} - 2} \right)\left( {1 - \chi} \right)}}.}$6. The method of claim 1 wherein the material azimuthal angle functionis a material Lode angle function.
 7. The method of claim 6 wherein thematerial Lode angle function is:$\mu_{0} = {\gamma + {\left( {1 - \gamma} \right){\left( \frac{6{\theta_{L}}}{\pi} \right)^{k}.}}}$8. The method of claim 1 wherein the material azimuthal angle functionis:$\mu_{0} = {{1 - {\left( {1 - \gamma} \right)\left( \frac{6{\theta }}{\pi} \right)^{k}\mspace{14mu} {for}\mspace{14mu} 0}} \leq \theta \leq {\frac{\pi}{6}.}}$9. The method of claim 1 wherein the material azimuthal angle functionis:$\mu_{0} = {{1 - {\left( {1 - \gamma} \right)\left( \frac{6{{{\pi/3} - \theta}}}{\pi} \right)^{k}\mspace{14mu} {for}\mspace{14mu} \frac{\pi}{6}}} < \theta \leq {\frac{\pi}{3}.}}$10. The method of claim 1 wherein the damage rule is a function of theratio of a plastic strain, ε_(p), to a fracture strain, ε_(f).
 11. Themethod of claim 10 wherein the damage rule is$D = {\left( \frac{ɛ_{p}}{ɛ_{f}} \right)^{m}.}$
 12. The method of claim10 wherein the damage rule is$D = {\frac{{\exp \left\lbrack {\lambda \left( \frac{ɛ_{p}}{ɛ_{f}} \right)}^{m} \right\rbrack} - 1}{{\exp \lbrack\lambda\rbrack} - 1}.}$13. The method of claim 1 wherein the damage rule is a function of theratio of a plastic distortion ε_(d) to a fracture strain ε_(p).
 14. Themethod of claim 13 wherein the damage rule is$D = {\left( \frac{ɛ_{d}}{ɛ_{p}} \right)^{m}.}$
 15. The method of claim13 wherein the damage rule is$D = {\frac{{\exp \left\lbrack {\lambda \left( \frac{ɛ_{p}}{ɛ_{f}} \right)}^{m} \right\rbrack} - 1}{{\exp \lbrack\lambda\rbrack} - 1}.}$16. The method of claim 1 wherein the weakening function is a linearfunction of the damage rule.
 17. The method of claim 1 wherein theweakening function is:w(D)=(1−D ^(β))^(1/η).
 18. The method of claim 17 wherein β is greaterthan one.
 19. The method of claim 17 wherein η is greater than one. 20.The method of claim 1 wherein an equivalent stress is used to determinematerial constitutive relationships of fracture, the equivalent stressbased on the weakening function and a matrix property of the ductilematerial.
 21. The method of claim 20 wherein the equivalent stress is:σ=w(D)σ_(M).
 22. A computer program product, tangibly embodied in aninformation carrier, the computer program product including instructionsbeing operable to cause a data processing apparatus to: determine amaterial pressure function describing dependence of a damage rule of atleast a portion of a ductile material on pressure; determine a materialazimuthal angle function describing dependence of the damage rule of theportion of the ductile material on principle stress; determine thedamage rule of the portion of the ductile material based on the materialpressure function and the material azimuthal angle function; determine aweakening function of the portion of the ductile material based on thedamage rule; and determine material constitutive relationships offracture in the portion of the ductile material based on the weakeningfunction and a matrix property of the ductile material.
 23. The computerprogram product of claim 22 wherein the instructions are incorporated inat least one of: a software or hardware package.
 24. The computerprogram product of claim 23 wherein the software package is LS-DYNA orABAQUS.
 25. A method of determining material constitutive relationshipsof damage accumulation in a ductile material undergoing plasticdeformation, the method comprising: determining a material pressurefunction describing dependence of a damage rule of at least a portion ofthe ductile material on pressure; determining a material azimuthal anglefunction describing dependence of the damage rule of the portion of theductile material on principle stress; and determining a nonlinear damageaccumulation rule for the portion of the ductile material based on thematerial pressure function and the material azimuthal angle function.26. The method of claim 25 further comprising: determining a weakeningfunction of the portion of the ductile material based on the nonlineardamage accumulation rule; and determining material constitutiverelationships of fracture in the portion of the ductile material basedon the weakening function and a matrix property of the ductile material.27. A method of determining damage accumulation characteristics ofmaterial constitutive relationships of fracture in a ductile materialundergoing plastic deformation, the method comprising: determining amaterial pressure function describing dependence of damage accumulationof at least a portion of the ductile material on pressure; determining amaterial azimuthal angle function describing dependence of damageaccumulation of the portion of the ductile material on principle stress;determining a damage rule of the portion of the ductile material basedon the material pressure function and the material azimuthal anglefunction; determining a damage accumulation in an integral form based onthe material damage rule, the pressure function, and the materialazimuthal angle function.
 28. The method of claim 27 further comprising:determining a weakening function of the portion of the ductile materialbased on the nonlinear damage accumulation rule; and determiningmaterial constitutive relationships of fracture in the portion of theductile material based on the weakening function and a matrix propertyof the ductile material.
 29. A method of determining materialconstitutive relationships of fracture in a ductile material undergoingplastic deformation, the method comprising: means for determining amaterial pressure function describing dependence of damage accumulationof at least a portion of the ductile material on pressure; means fordetermining a material azimuthal angle function describing dependence ofdamage accumulation of the portion of the ductile material on principlestress; means for determining a damage rule of the portion of theductile material based on the material pressure function and thematerial azimuthal angle function; means for determining a weakeningfunction of the portion of the ductile material based on the damagerule; and means for determining material constitutive relationships offracture in the portion of the ductile material based on the weakeningfunction and a matrix property of the ductile material.
 30. A method ofdetermining material constitutive relationships of fracture in a ductilematerial undergoing plastic deformation, the method comprising:determining a material pressure function describing dependence of damageaccumulation of at least a portion of the ductile material on pressure;determining a material azimuthal angle function describing dependence ofdamage accumulation of the portion of the ductile material on principlestress; determining a damage rule of the portion of the ductile materialbased on the material pressure function and the material azimuthal anglefunction; determining a weakening function of the portion of the ductilematerial based on the damage rule, wherein the weakening function is afunction of damage accumulation and the rate of change of the weakeningfunction with respect to damage accumulation is dependent on damageaccumulation; and determining material constitutive relationships offracture in the portion of the ductile material based on the weakeningfunction and a matrix property of the ductile material.